US  Army  Corps 
of  Engineers 


WAVE  INFORMATION  STUDIES 
OF  US  COASTLINES 

WIS  REPORT  27 


USER’S  GUIDE  TO  THE  WAVE  INFORMATION 
STUDIES  (WIS)  WAVE  MODEL,  VERSION  2.0 


by 

J.  M.  Hubertz 

Coastal  Engineering  Research  Center 

DEPARTMENT  OF  THE  ARMY 
Watenways  Experiment  Station,  Corps  of  Engineers 
3909  Halls  Ferry  Road,  Vicksburg,  Mississippi  39180-6199 


rx) 


June  1992 
Final  Report 


Approved  For  Public  Release;  Distribution  Is  Unlimited 


Prepared  for  DEPARTMENT  OF  THE  ARMY 
US  Army  Corps  of  Engineers 
Washington,  DC  20314-1000 


Destroy  this  report  when  no  longer  needed.  Do  not  return 
it  to  the  originator. 


The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position  unless  so  designated 
by  other  authorized  documents. 


The  contents  of  this  report  are  not  to  be  used  for 
advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of 
such  commercial  products. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  NO.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  t  hour  per  response,  including  the  time  for  reviewing  instructions.  searchir>g  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  cottection  of  information.  Send  comments  regarding  this  burden  estimate  Of  any  other  aspect  of  this 
colleaion  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway.  Suite  1204.  Arlington.  VA  22202'4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0704*0188),  Washir>gton,  DC  20503. 


1.  AGENCY  USE  ONLY  (leave  b/anfc) 


4.  TITLE  AND  SUBTITLE 


REPORT  DATE 

June  1992 


3.  REPORT  TYPE  AND  DATES  COVERED 

Final  renort 


5.  FUNDING  NUMBERS 


User’s  Guide  to  the  Wave  Information 
Studies  (WIS)  Wave  Model,  Version  2.0 


6.  AUTHOR(S) 

J.  M.  Hubertz 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

USAE  Waterways  Experiment  Station,  Coastal  Engineering 
Research  Center,  3909  Halls  Ferry  Road,  Vicksburg,  MS 
39180-6199 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


WIS  Report  27 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

US  Army  Corps  of  Engineers 
Washington,  DC  20314-1000 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 


Available  from  National  Technical  Information  Service,  5285  Port  Royal  Road 
Springfield,  VA  22161. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution  is  unlimited. 


13.  ABSTRACT  (Maximum  200  words) 

This  report  provides  guidance  on  application  of  the  Wave  Information  Studies  Wave  Model, 
Version  2.0  (WISWAVE  2.0),  and  a  description  of  the  upgrades  from  Version  1.0.  The  present 
version  operates  in  all  water  depths,  while  the  previous  version  operated  only  in  depths  sufficient  not 
to  affect  wave  propagation.  The  present  version  also  allows  for  changing  water  depths  during  a 
simulation  as  would  be  the  case  in  storm  surge  applications  or  regions  of  large  tidal  amplitude. 

The  structure  and  operation  of  the  model  are  described.  Input  necessary  to  operate  the  model 
is  described  and  an  example  input  data  set  is  provided.  The  theory,  equations,  and  solution 
techniques  contained  within  the  model  are  described  in  the  references  provided.  The  theory, 
equations,  and  solution  techniques  associated  with  the  upgrade  to  include  shallow-water  effects  are 
described  in  Appendix  A.  Documented  subroutines  associated  with  the  upgrade  are  provided  in 
Appendix  B. 

(Continued) 


14.  SUBJECT  TERMS 

Numerical  model 
User’s  guide 
Water  waves 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

UNCLASSMED  UNCLASSIFIED 


15.  NUMBER  OF  PAGES 

42 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


Standard  Form  298  (Rev  2-89) 

P^e^crib^d  by  anSi  Std  739->8 
298-102 


13.  (Concluded). 


Comparisons  of  model  results  to  measurements  are  provided  in  WIS  reports  on  the  Great 
Lakes  (Reports  22-26)  and  the  Atlantic  hindcast  (WIS  Report  30).  WISWAVE  1.0  was  used  to 
produce  the  Great  Lakes  wave  hindcasts  described  in  WIS  Reports  22-26.  WISWAVE  2.0  was  used 
to  produce  the  Atlantic  hindcast,  WIS  Report  30. 


PREFACE 


The  Wave  Information  Studies  (WIS)  were  authorized  by  Headquarters,  US 
Army  Corps  of  Engineers  (HQUSACE)  in  1976  at  the  US  Army  Engineer  Waterways 
Experiment  Station  (WES)  to  produce  a  wave  climatology  for  US  coastal  waters. 
The  WIS  is  part  of  the  Coastal  Field  Data  Collection  Program  managed  by  the 
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Technical  Monitors  are  Messrs.  John  H.  Lockhart,  J.  Housley,  F.  Campbell,  and 
J .  Crews . 

This  report,  the  27th  in  a  series,  provides  guidance  on  the  operation  of 
the  latest  version  of  the  WIS  wave  model  (WISWAVE  2.0)  developed  under 
contract  to  Dr.  Donald  T.  Resio  of  Offshore  and  Coastal  Technology 
Incorporated  (OCTI) .  The  structure  and  operation  of  the  model  are  described 
and  example  input  is  provided.  The  following  members  of  the  WIS  group  in  the 
Coastal  Oceanography  Branch  (COB),  CERC,  were  instrumental  in  documenting  and 
verifying  the  computer  code:  Dr.  J.  M.  Hubertz,  Mr.  D.  B.  Driver,  Mses.  B.  A. 
Tracy,  R.  D.  Reinhard,  and  W.  A.  Brandon.  Dr.  Resio  provided  valuable 
assistance  during  documentation  and  testing  of  the  model. 

The  work  was  conducted  under  the  direct  supervision  of  Dr.  Martin  C. 
Miller,  Chief,  COB,  and  Mr.  H.  L.  Butler,  Chief,  Research  Division,  CERC;  and 
under  general  supervision  of  Dr.  J.  R.  Houston  and  Mr.  C.  C.  Calhoun,  Jr., 
Director  and  Assistant  Director,  CERC,  respectively.  The  word  processing  of 
this  report  was  done  by  Ms .  M.  J.  Stauble,  COB.  The  report  was  edited  by  Ms. 
Janean  C.  Shirley,  Information  Technology  Laboratory,  WES. 

At  the  time  of  publication  of  this  report.  Director  of  WES  was 
Dr.  Robert  W.  Whalin.  Commander  and  Deputy  Director  was  COL  Leonard  J 
Hassell,  EN. 
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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC)  UNITS  OF  MEASUREMENT 


Non- SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI 
(metric)  units  as  follows: 

_ Multiply _  _ _  To  Obtain 

degrees  (angle)  0.01745329  radians 
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USER'S  GUIDE  TO  THE  WAVE  INFOEmATION  STUDIES  (WIS)  WAVE  MODEL 


VERSION  2,0 

PART  I :  INTRODUCTION 

1.  A  number  of  numerical  wave  models  have  been  used  under  the  Wave 
Information  Studies  (WIS)  to  produce  hindcast  wave  information  along  US 
coastlines.  The  earliest  version  (Resio  and  Vincent  1977,  1978),  was  used  to 
produce  hindcasts  of  storms  on  the  Great  Lakes.  The  next  version  (Resio  1981) 
was  used  for  the  Atlantic  Coast  Phases  I  and  II  hindcasts.  The  Pacific  Coast 
hindcasts.  Phases  I  and  II,  were  produced  by  a  version  described  by  Resio  and 
Tracy  (1983)  and  the  Gulf  of  Mexico,  Phase  II  hindcast  by  a  version  described 
by  Resio  (1987  and  1988). 

2.  Each  successive  model  improved  in  some  way  over  its  predecessor. 
However,  each  model  was  more  a  separate  entity  than  an  improvement  to  a  basic 
model  framework.  The  present  version  (Resio  and  Perrie  1989)  is  structured  to 
allow  improvements  to  a  basic  framework  so  that  users  see  little  change  in 
applying  the  model  as  improved  versions  are  developed. 

3.  The  first  version  of  the  present  model  (hereafter  referred  to  as 
WISWAVE  1.0)  was  used  to  produce  the  Great  Lakes  hindcasts  and  is  described  in 
each  of  the  WIS  reports  on  the  Great  Lakes  (e.g.,  Hubertz,  Driver,  and 
Reinhard  1991).  The  present  version  (hereafter  referred  to  as  WISWAVE  2.0) 
operates  in  all  water  depths,  while  WISWAVE  1.0  only  operated  in  depths  deep 
enough  not  to  affect  wave  propagation.  The  WISWAVE  2.0  also  allows  for 
changing  water  depths  during  a  simulation  as  would  be  the  case  in  a  storm 
surge,  application  or  regions  of  large  tidal  amplitude. 

4.  The  WISWAVE  2.0  was  used  to  re-hindcast  wave  conditions  along  the 
Atlantic  Coast  (Hubertz  et  al.,  in  preparation)  for  the  period  1956-75  to 
improve  previous  hindcast  results.  The  winds  developed  for  the  original  WIS 
hindcasts  (Corson,  Resio,  and  Vincent  1980;  Resio,  Vincent,  and  Corson  1982) 
were  used  in  the  re-hindcast  since  they  provide  an  accurate  representation  of 
the  wind  climatology.  Present  computer  technology  allowed  greater  resolution 
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(l.O-deg*  and  0.25-deg  spacing  for  Phases  I  and  II,  respectively,  versus  the 
2.0-deg  and  0.5-deg  spacing  used  in  the  original  hindcasts) . 

5.  The  WISWAVE  model  is  written  using  standard  FORTRAN  commands.  It 
can  be  executed  on  computers  ranging  from  microcomputers  to  super  computers. 
However,  the  model  employs  multi-dimensioned  arrays  dependent  on  the  number  of 
grid  points  covering  a  geographical  area,  and  the  number  of  frequency  and 
angle  bands  specified  to  describe  the  directional  wave  spectra.  Thus,  super 
computers  are  needed  to  perform  simulations  on  large  bodies  of  water  due  to 
the  large  number  of  calculations  involved.  The  version  of  WISWAVE  that  WIS 
uses  has  been  optimized  for  execution  on  the  CRAY  YM-P  super  computer  at  the 
US  Army  Engineer  Waterways  Experiment  Station.  The  structure  and  flow  control 
of  the  model  are  discussed  next. 


*i\  table  of  factors  for  converting  non-Sl  units  of  measurement  to 
SI  (metric)  unit.s  is  presented  on  page  3. 
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PART  II:  STRUCTURE  AND  OPERATION 


6.  There  are  two  main  DO  loops  in  the  model.  One  (the  inner)  is 
contained  within  the  other  (the  outer).  The  inner  loop  steps  through  time 
steps  from  the  beginning  of  a  simulation  to  the  end  and  calculates  wave 
conditions  at  each  point  on  a  latitude -longitude  grid  as  a  function  of  input 
conditions.  The  outer  loop  reads  input  as  a  function  of  time.  This  input  can 
consist  of  winds  over  the  model  grid,  wave  conditions  on  the  grid  boundary,  or 
water  level  information  over  the  grid.  This  information  is  specified  at  time 
intervals  during  the  simulation  to  update  changing  wind,  wave,  and  water  level 
conditions . 

7.  Prior  to  entering  these  loops,  a  number  of  subroutines  read  required 
options  for  the  particular  application  and  initialize  values  of  constants  and 
variables.  Information  is  output  as  simulated  time  progresses  and  information 
is  saved  if  subsequent  runs  of  the  same  application  are  required.  A  simple 
flow  diagram  is  shown  in  Figure  1. 

Description  of  Subroutines 

8.  Subroutine  RDOPT  (Read  Options)  reads  the  model  input  options 
describing  an  application  from  unit  20.  This  includes  a  land/water  matrix 
which  identifies  each  point  in  the  latitude-longitude  grid  as  land  (0)  or 
water  (1),  and  identifies  boundaries  of  a  nested  grid  if  used.  A  nested  grid 
is  a  subarea  within  the  original  grid  boundary.  It  is  defined  by  grid  lines 
of  the  first  level  grid,  but  has  a  higher  density  of  grid  lines,  which  allows 
greater  resolution  of  wave  information  in  space .  The  initial  depths  at  each 
grid  point,  if  depths  are  used,  are  also  read  from  this  file.  Updates  to  the 
depth  field  as  a  function  of  time  are  read  from  unit  14  when  wind  fields  are 
updated . 

9.  Subroutine  CKSTR  (Check  Start)  either  reads  (from  unit  23)  or  writes 
(to  unit  16)  information  to  be  used  in  a  warm  start.  A  warm  start  is  a  run 
which  uses  results  from  the  end  of  a  previous  run  to  initialize  variables  in 
the  model.  This  allows  computations  to  proceed  from  where  they  left  off 
without  a  "spin-up."  If  a  warm  start  is  not  used,  there  will  be  some  time 
involved  for  wave  conditions  to  develop  for  the  given  winds.  The  length  of 
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WISWAV 


END 


Figure  1.  Flow  diagram  of  WISWAVE  2.0 

time  is  related  to  the  size  of  the  water  body,  and  can  be  on  the  order  of  a 
few  days  for  an  ocean. 

10.  Subrovitine  I"'",.-  (Initialize  Boundary  Information)  is  only  used 
when  the  present  run  is  using  boundary  information  saved  from  a  previous  run. 
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That  is,  a  grid  nested  within  a  previously  used  larger  grid  is  being  used.  It 
checks  the  land/water  matrix  to  identify  the  location  of  points  where  boundary 
information  was  saved  and  interpolates  and  extrapolates  this  information  to 
the  grid  points  along  the  boundary  of  the  present  grid  being  used. 

11.  Subroutine  RDBND  (Read  Boundary  Information)  is  only  used  when  the 
present  run  is  using  boundary  information  saved  from  a  previous  run.  It  reads 
(from  unit  25)  the  previously  saved  boundary  information  (written  by 
subroutine  WBOUND  on  unit  17)  and,  using  the  weighting  factors  calculated  in 
INITB ,  calculates  boundary  conditions  along  the  model  input  boundary. 

12.  Subroutine  RDWIN  (,?.ead  Wind  Input)  is  used  to  read  wind  speeds  and 
directions  (from  unit  21)  for  all  grid  points  at  fixed  time  intervals  during 
the  simulation.  Wind  speeds  and  directions  are  interpolated  to  the  timestep 
interval,  from  times  when  wind  fields  are  supplied.  The  time  step  is  the 
interval  in  time  between  which  wave  conditions  are  calculated  over  the  grid. 

It  must  be  less  than  the  time  needed  for  wave  energy  to  travel  between 
adjacent  grid  points. 

13.  Subroutine  INITSPR  (Initialization  for  Shallow-Water  Propagation) 
is  used  to  initialize  all  parameters  dependent  on  depth.  It  is  called  once 
for  a  shallow- water  application.  It  is  also  called  every  time  the  depths  are 
updated  as  in  a  hurricane  surge  or  application  with  a  large  tide  range. 

14.  Subroutine  DPROP  (Deep-Water  Propagation)  is  used  to  propagate  wave 
energy  over  the  grid  assuming  deep  water.  Water  depths  are  set  at  999  m 
everywhere,  and  refraction  and  shoaling  do  not  take  place. 

15.  Subroutine  SPROP  (Shallow-Water  Propagation)  is  used  to  propagate 
wave  energy  over  the  grid  accounting  for  the  effects  of  shallow  water 
(refraction,  shoaling). 

16.  Subroutine  FLXSTD  (Flux  of  Energy  through  the  Spectrum)  is  the 
subroutine  that  calculates  the  two-dimensional  spectrum  at  each  grid  point 
based  on  the  wind  input  and  transfer  of  energy  within  the  spectrum  due  to  the 
wave  interaction  source  term. 

17.  Subroutine  WBOUND  (Write  Boundary  Information)  writes  information 
to  a  file  (on  unit  17)  for  later  use  in  a  nested  grid  application.  This  is 
the  same  file  subroutine  RDBND  accesses  (on  unit  25)  when  it  is  reading  in 
lnf(  rmation  on  the  boundary  of  a  nested  grid. 
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18.  Subroutine  NDATE  (Calculate  Date)  calculates  the  date  used  on 
output  information  from  the  beginning  date  and  length  and  number  of  time  steps 
between  specified  output  times. 

19.  Subroutine  OUTP  (Write  Output  Information)  writes  out  the  results 
of  the  calculations  (on  unit  93)  at  the  times  specified  in  the  options  file  by 
the  parameter  NTMS . 


Model  Input 


20.  The  grid  system  used  by  the  model  is  the  system  of  latitude  and 
longitude  lines  on  the  Earth's  surface.  Lines  of  latitude  are  denoted  by 
rows,  referenced  by  values  of  J,  and  lines  of  longitude  are  columns, 
referenced  by  values  of  I.  The  origin  of  a  grid  is  in  the  lower  left-hand 
corner  (Figure  2).  For  small  geographical  areas  on  the  Earth's  surface,  such 
as  one  of  the  Great  Lakes,  the  distance  between  specified  grid  lines  can  be 
approximated  as  equal  over  the  grid.  For  larger  areas,  the  curvature  of  the 
Earth  is  taken  into  account  to  calculate  the  distance  between  grid  lines  so 
that  waves  are  propagated  over  the  proper  distances. 


LONGITUDE 


COLUMNS  I  =  1,  NX 

Figure  2.  Schematic  of  model  grid 


21.  The  dimensions  of  the  arrays  in  the  model  are  specified  by  two 
parameter  files  PARAl . INC  and  PARA2 . INC  which  are  included  in  the  program 
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during  compilation.  Thus,  they  must  be  prepared  and  available  in  the  same 
directory  as  the  source  code  before  the  program  is  compiled.  Examples  of 
these  files  are  shown  below  for  the  sample  case  discussed  later. 

PARAl.INC  parameter  ( idmn-=18 , j dmn=35 , if=20 , ia=16 , nobpts=iO , nbn=30) 

PARA2 . INC  parameter  ( idmn=18 , j dran=35 , if=20 , ig=16 , nobpts=10 , nbn=30) 

The  dimensions  for  the  grid  size,  columns  and  rows,  are  set  with  (idmn)  and 
(jdmn)  respectively.  The  dimensions  for  frequency  and  direction  increments 
are  set  with  parameters  (if)  for  frequency  and  both  (ia)  and  (ig)  for 
direction.  If  boundary  information  is  saved,  (nobpts)  sets  the  dimensions  for 
the  number  of  save  boundary  points.  If  an  application  is  making  use  of 
previously  saved  boundary  Information,  (nbn)  specifies  the  number  of  boundary 
points  in  the  present  application  which  will  be  made  up  of  interpolated  and 
extrapolated  points  determined  from  those  previously  saved. 

22.  There  are  a  number  of  parameters  that  specify  the  details  of  a 
particular  application.  These  are  given  in  an  options  file  which  is  made  up 
of  a  number  of  lines  or  records.  Values  of  parameters  need  to  be  specified  in 
the  order  shown  below  for  each  record  and  separated  by  at  least  one  space. 

Note  that  some  parameters  are  integers  while  others  are  real  numbers.  A 
summary  of  input  and  output  unit  numbers  and  associated  information  is 
provided  in  Table  1. 


LOGICAL  RECORD  1 


NX  NY  NANG  NFRQ  DL  DT  MSTA  INPLEV 


NX 

NY 

NANG 

NFRQ 

DL 

DT 

MSTA 

INPLEV 


The  number  of  columns  in  the  grid. 

The  number  of  rows  in  the  grid. 

The  number  of  angle  bands  in  the  discrete  spectrum. 

The  number  of  frequency  bands  in  the  discrete  spectrum. 
The  distance  between  grid  points  in  kilometers. 

The  time  step  in  seconds. 

The  number  of  points  at  which  results  should  be  saved. 
Option  to  specify  input  level  of  winds. 

0  =  10  m  above  sea  surface. 

1  =  20  m  above  sea  surface. 
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Table  1 

Input  and  Output  Unit  Numbers  and  Associated  Information 


Unit  Number 


20 

23 

25 

21 


Unit  Number 


93 

16 

17 


Input  Files 

_ Data 


Options  file  (see  page  8) 

Warm  start  information  (see  page  8) 

Boundary  data  for  nested  grid  (see  page  10) 
Wind  fields  (see  p.  10) 

Output  Files _ 

_ Data _ 


One  record  summary  at  output  times 
(see  page  11  for  description) 

Warm  start  information  (see  page  8) 

Nested  grid  boundary  information 
(see  page  10) 


LOGICAL  RECORD  2 

NSTR  NORD  NTMS  NHR  MXHR  IBND  ISTDEP 


NSTR 


NORD 


Option  to  specify  whether  results  from  a  previous  run 
will  be  used  to  start  present  run  (warm  start) . 

Cold  start,  no  restart  data  read. 

Cold  start,  but  save  results  for  a  warm  start. 
Warm  start,  but  do  not  save  results  for 
another  warm  start. 

Warm  start,  and  save  results  for  another  warm 
start . 

read  in  windfields  or  not. 

Read  in  windfields. 

Do  not  read  in  windfields,  but  specify  fixed 
speeds  and  directions. 


0 
1 
2 

3  = 

Option  to 
0  = 
1  = 
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NTMS  Number  of  time  steps  between  output  of  results.  For 

example,  output  Is  desired  every  3  hr  and  time  step  is  30  min  (0.5 
hr) ,  NTMS  =  3/0.5  =  6 

NHR  Number  of  hours  between  input  of  winds . 

MXHR  Maximum  number  of  wind  inputs  in  run.  This  value 

should  be  consistent  with  information  in  record  6. 

IBND  Option  to  read  and  write  results  on  boundary  of  grid 

for  use  in  nested  grid  applications. 

0  =  No  boundary  data  read  or  written. 

1  =  Boundary  data  written,  but  none  read. 

2  =  Boundary  data  read,  but  none  written. 

3  =  Boundary  data  read  and  written. 

ISTDEP  Option  to  use  depths  or  assume  uniform  value. 

0  =  Depths  input  (meters) 

1  =  Depth  assumed  constant  (999  m) . 

LOGICAL  RECORD  3 
DLAT  XLATO  ICURV  lUPDAT 

DLAT  Number  of  degrees  of  latitude  between  rows  of  grid. 

XLATO  Southernmost  latitude  (degrees)  of  grid.  Negative  if  in 

Southern  Hemisphere. 

ICURV  Option  to  include  curvature  effects  of  Earth  in 

propagation. 

0  =  Do  not  include. 

1  =  Include. 

IDPRP  Option  to  choose  deep  or  shallow-water  propagation. 

0  =  Shallow  water 
1  =  Deep  water 

lUPDAT  Option  to  update  still -water  level  during  simulation. 

0  =  Do  not  update. 

1  =  Update  water  levels  at  frequency  of  wind 
input , 

LOGICAL  RECORD  h 

FREQ(K)  Mid-band  values  of  frequency  for  NFRQ  bands.  Can  be  on 

more  than  one  line. 

LOGICAL  RECORD  5  (Only  needed  if  MSTA  >  0) 
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K  lOUT(K) ,  JOUT(K) 


K  Sequential  number  of  output  point  (station  number) 

lOUT(K)  Grid  column  at  station  K. 

JOUT(K)  Grid  row  at  station  K. 

LOGICAL  RECORD  6 

IDFIRST  IDLAST  IDBNDl  IDBND2 

IDFIRST  The  date/time  code  to  start  a  simulation  (for  example, 

90010100  signifies  1990  January  1,  00  hr). 

IDLAST  The  date/time  code  to  end  a  simulation. 

IDBNDl  The  date/time  code  to  begin  saving  boundary 

information. 

IDBND2  The  date/time  code  to  end  saving  boundary  information. 

LOGICAL  RECORD  7  (Only  needed  if  ISTDEP=0) 

DEP(I,J)  This  array  contains  the  depth  at  each  grid  point  if  ISTDEP=0. 

Values  should  be  positive,  in  meters,  and  arranged  for  all  columns 
starting  with  row  NY  (see  example  in  sample  application) .  Format 
statement  401  in  RDOPT  can  be  adjusted  to  meet  the  user's  needs. 

LOGICAL  RECORD  8 

IB0UND(I,J)  A  matrix  of  values  defining  whether  a  grid  point  is  land,  water, 
or  a  boundary  point  for  a  nested  grid  application.  Values  are 
arranged  for  input  in  the  same  manner  as  the  depths.  The 
following  numerical  values  are  used  to  designate  each  grid  point. 
Examples  are  shown  in  the  sample  application.  Format  is  8011. 

0  -  Land  Point 

1  -  Water  Point 

2  -  Extrapolated  Boundary  Point  (at  beginning  of  boundary 

line) 

3  -  Input  Boundary  Point  (coincident  point  in  large  and 

nested  grid) 

4  -  Interpolated  Boundary  Point  (interpolated  at  points  on 

boundary  of  nested  grid  between  (3)  type  points) 

5  -  Extrapolated  Points  (at  end  of  boundary  line) 

LOGICAL  RECORD  9  Only  needed  if  IBND  is  1  or  3 . 

NBPO  The  number  of  boundary  points  at  which  information  is 

output  for  a  nested  grid  run. 
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LOGICAL  RECORD  10 


(Only  needed  if  IBND  is  1  or  3) 

IBPO(K) ,  K— l.NBPO  Ordered  I  locations  of  boundary  output  points 
from  beginning  of  boundary  curve  to  end. 

LOGICAL  RECORD  11  (Only  needed  if  IBND  is  1  or  3) 

JBPO(K) ,  K”1,NBP0  Ordered  J  locations  of  boundary  output  points 
from  beginning  of  boundary  curve  to  end. 

LOGICAL  RECORD  12  (Only  needed  if  IBND  is  2  or  3) 

ILS,  JLS  The  I,J  location  (in  the  nested  grid  coordinates)  of  the  boundary 
point  at  the  beginning  of  the  offshore  boundary  curve  in  the 
nested  grid.  It  will  usually  be  an  IBOUND(I , J)=2  type  point. 

LOGICAL  RECORD  13  (Only  needed  if  NORD=l) 

lUWS,  lUDIR,  ISHFT,  NSHFT,  IWNDN,  IDEGN 


lUWS 

lUDIR 

ISHFT 

NSHFT 

IWNDN 


A  constant  value  of  wind  speed  (m/sec)  over  the  whole 
grid. 

The  wind  direction  in  degrees.  The  wind  angle  is  defined  in 
the  compass  sense,  i.e.,  from  North  =  0  deg,  from  East  =  90 
deg,  etc. 

Option  to  change  wind  speed  and  direction  in  test  case. 
Number  of  hours  in  simulation  prior  to  shift  in  wind. 

Wind  speed  after  shift. 


IDEGN  Wind  direction  after  shift. 


23.  Two  files  need  to  be  prepared  by  the  user.  One  is  the  options 
file,  and  the  other  the  wind  input  file.  The  parameters  in  the  options  file 
are  described  above.  The  wind  speed  and  direction  at  every  grid  point  in  the 
grid  are  read  in  every  NHR  hours  in  the  simulation.  The  wind  file  is  composed 
of  a  series  of  wind  fields,  each  consisting  of  a  header  record  containing  the 
date/time  code  for  the  following  winds.  Next  the  wind  speeds  (meters/second) 
are  read,  in  the  same  order  as  the  depths  and  land/water  boundary  matrix,  and 
then  the  wind  directions,  also  in  the  same  order. 

24.  Wind  directions  should  be  in  the  meteorological  or  compass  sense, 
i.e.,  0  deg  is  a  wind  blowing  from  the  north,  90  deg,  from  the  east,  etc. 
Winds  are  read  in  by  subroutine  RDWIN.  The  present  formats  for  the  header 
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record  and  speed  and  direction  are  I  10,  F5.1,  and  F5.0,  respectively.  The 
user  may  adjust  the  format  and  make  speed  and  direction  conversions  as  needed 
in  this  subroutine,  or  reformat  existing  wind  fields  to  conform  to  the  above 
format  and  conventions. 

25.  Wind  speeds  (m/sec)  should  be  representative  of  either  a  10-  or 
20 -m  elevation  above  the  sea  surface.  The  particular  level  is  specified  by 
the  parameter  INPLEV  in  record  1  of  the  options  file.  If  the  winds  are  not  at 
one  of  these  levels,  but  at  some  other  level  between  0  and  20  m,  they  can  be 
adjusted  to  one  of  these  levels  using  Equation  3-26,  page  3-26  of  the  Shore 
Protection  Manual  (SPM)  (US  Army  Engineer  Waterways  Experiment  Station  1984) . 

Model  Output 

26.  The  model  calculates  the  discretized  two-dimensional  (frequency  by 
direction)  wave  spectrum  at  every  time  step  at  every  grid  point  during  the 
simulation.  Potentially,  all  of  this  information  could  be  output,  but  would 
be  overwhelming  to  comprehend.  Generally  a  user  is  interested  in  the  basic 
wave  parameters  as  a  function  of  time  at  particular  points.  These  parameters 
are:  (a)  wave  height,  (b)  peak  or  mean  wave  period,  and  (c)  mean  wave 
direction.  These  values,  along  with  a  station  number  (location  identifier), 
the  date/time  code,  the  wind  speed  and  direction,  and  water  level  with  respect 
to  a  datum  are  output  on  unit  93. 

27.  Other  information  which  describes  the  wave  conditions  in  more 
detail  is  also  provided  on  this  output  record.  This  consists  of  the  percent 
of  total  energy  in  each  frequency  band  and  the  mean  direction  of  each 
frequency  band.  Additional  space  is  also  available  on  this  record  for  the 
speed  and  direction  of  flow  of  currents  at  the  chosen  output  location.  This 
is  in  anticipation  of  a  future  upgrade  to  add  currents  to  the  wave  model. 

28.  An  example  of  the  output  record  is  shown  below. 

OOl  i>0102S0f.  IK4  11  10  165  06  16.5  00  00  00  001*46  0  0  0  0  0  0  1  3  61311  7  5  8IOI3I2  63  lOOOOOO  04646*646  0  0  0  0''0000 

The  first  three  digits  are  the  station  number.  The  next  eight  are  the  date/ 
time  group,  which,  taken  in  four  groups  of  two  represent  the  year,  month,  day, 
and  hour.  The  example  above  represents  0600  on  October  28,  1990.  The  next 
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three  digits  represent  the  significant  wave  height  in  meters  times  ten.  In 
this  example,  the  value  is  2.4  m.  The  next  two  values  are  the  peak  and  mean 
wave  periods  to  the  closest  second.  In  this  example,  they  are  13  and  10  sec, 
respectively.  The  next  three  digits  represent  the  mean  wave  direction  in 
degrees  from  which  the  waves  are  coming  with  respect  to  compass  directions. 
That  is,  0  deg  is  from  the  north,  90  deg,  from  the  east,  etc.  In  this 
example,  the  waves  are  coming  from  165  deg  or  from  a  direction  between  SSE  and 
S.  The  next  two  digits  are  the  wind  speed  in  m/sec,  6  m/sec  in  this  case. 

The  next  three  are  the  wind  direction  in  degrees  in  the  same  convention  as  the 
waves.  In  this  example,  both  the  wind  and  waves  are  coming  from  165  deg.  The 
next  two  digits  represent  the  water  level  above  a  datum.  There  will  be  a  non¬ 
zero  value  here  only  if  the  water  level  is  changed  during  the  simulation  as  in 
a  hurricane  surge  or  large  tidal  range  case.  The  value  is  in  meters  with 
respect  to  the  datum  chosen;  for  example,  mean  sea  level  or  mean  low  water,  or 
some  other  datum.  The  next  two  digits  represent  the  speed  of  current  flow 
along  a  fixed  latitude  or  in  the  x  direction  (see  Figure  2).  Values  are  in 
meters/second  times  10.  The  next  two  digits  are  the  similar  quantity  for  flow 
in  the  y  direction.  Positive  values  are  toward  the  east  and  north, 
respectively.  The  effect  of  currents  on  wave  propagation  and  source  terms  is 
not  accounted  for  in  WISWAVE  2.0,  but  is  planned  for  the  future.  The  next  six 
digits  represent  the  total  energy  in  the  discrete  part  of  the  spectrum 
multiplied  by  10,000  and  in  units  of  meters  squared. 

29.  The  next  40  spaces  contain  the  percent  of  total  energy  in  the 
discrete  part  of  the  spectrum  contained  in  each  of  the  20  frequency  bands  for 
all  directions.  The  center  frequencies  associated  with  each  of  these  bands 
are  given  in  record  4  of  the  input  options  file.  The  hindcast  frequency 
spectrum  is  obtained  using  these  frequencies,  the  total  discrete  energy,  and 
the  percentages  in  each  band. 

30.  The  last  40  spaces  contain  the  mean  directions  in  each  frequency 
band,  expressed  as  a  percent  of  360  deg.  In  the  example  above,  in  frequency 
band  eight,  the  46  represents  46  percent  of  360  deg  or  approximately  165  deg. 
Thus,  in  this  frequency  band,  the  merin  direction  from  which  waves  are  coming 
is  about  165  deg  or  from  the  SSE.  The  mean  across  all  frequency  bands  should 
agree  with  the  mean  direction  of  the  spectrum  given  by  the  sixth  group  of 
numbers  in  the  record.  There  may  be  some  difference  due  to  rounding  off  in 
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the  averaging  process  and  the  resolution  of  4  deg  imposed  by  expressing  the 
directions  as  a  percentage  of  360  deg. 

31.  Additional  output  may  be  obtained  by  referring  to  subroutine  OUTP 
and  modifying  the  coding  to  obtain  what  is  desired.  Examples  are  the  complete 
two-dimensional  frequency  by  direction  spectrum,  or  values  of  wave  parameters 
such  as  height,  period,  and  direction  over  the  whole  grid  at  particular  times. 

Sample  Input  Files 

32.  The  following  sample  input  files  are  provided  for  the  user  to  he_p 

in  understanding  input  to  the  model.  The  example  is  for  Lake  Michigan.  The 

numerical  model  grid  for  this  application  is  shown  in  Figure  3.  An 

explanation  of  each  of  the  input  parameters  is  provided  below  for  this 

application.  A  general  description  is  provided  above  in  the  section  titled 

"Model  Input." 
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1  0  6  3  40  1  0 
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33.  The  first  two 
numbers  in  the  first  line  or 
record  of  input  are  18  and  35. 
For  this  application,  these 
numbers  indicate  that  there  are 
18  columns  (NX=18)  or  values  of 
the  I  index  and  35  rows  (NY=35) 
or  values  of  the  J  index  (see 


Figure  3)  . 


There  are  16 


direction  bands  (NANG=16)  ,  each 


of  which  is  22.5  deg  in  width. 
These  correspond  to  the  points 
of  the  compass,  i.e.,  N,  NNE, 
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Figure  3.  Example  of  model  grid 
of  the  compass,  i.e.,  N,  NNE,  for  fake  Michigan 

NE,  ENE,E,  etc. ,  and  are  typically  the  number  of  direction  bands  specified  in  an 
application.  There  are  20  frequency  bands  specified  (NFRQ=20) ,  which  is  a 
typical  number  to  use.  The  spacing  between  grid  lines  (DL)  is  16.1  km,  which 
corresponds  to  10  miles.  The  time  step  (DT)  is  1,800  sec  or  30  min.  Information 
at  one  grid  intersection  point  is  being  output  (MSTA=1) .  The  winds  being  input 
are  applicable  to  an  elevation  10  m  above  the  sea  surface  (INPLEV=0) . 
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34.  The  first  value  in  the  second  line  or  record  of  input  indicates 

that  this  is  the  initial  run  for  this  application  (NSTR=1) ;  that  is,  no 
previous  results  are  available  as  a  warm  start,  but  information  for  a  warm 
start  will  be  saved  at  the  end  of  this  run.  Wind  speeds  and  directions  will 
be  read  in  for  each  grid  point  for  successive  times  (NORD=0) .  Results  will  be 

output  every  6  time  steps  (NTMS=6)  or  in  this  case  every  3  hr,  since  a  time 

step  equals  0.5  hr.  New  wind  information  will  be  input  every  3  hr  (NHR=3)  and 
a  total  of  40  wind  fields  (MXHR=40)  will  be  input.  This  corresponds  to  120  hr 

of  simulation,  or  5  days.  Boundary  information  will  be  written  at  a  number  of 

grid  points  defining  the  outer  boundary  of  a  grid  nested  within  this  grid 
(IBND=1).  Depths  will  be  input  (ISTDEP=0)  to  account  for  refraction, 
shoaling,  and  depth  limitation  of  wave  energy. 

35 .  The  third  record  begins  with  the  increment  in  degrees  of  latitude 
(DLAT)  between  rows  of  the  grid,  whose  value  in  this  case  is  0.145  deg.  The 
southern-most  latitude  of  the  grid  (XLATO)  is  41.55  deg.  The  size  of  the  grid 
on  the  Earth's  surface  is  relatively  small,  so  curvature  effects  of  the 
Earth's  surface  are  not  included  (ICURV=0).  Shallow  water  is  present,  so 
IDPRP=0.  Water  levels  are  not  being  changed  during  the  simulation,  so 
IUPD^T-0) . 

36.  The  fourth  record  contains  the  NFRQ  values  of  frequency  of  the 
discretized  frequency  spectrum.  These  are  the  values  at  the  center  of  the 
bands.  The  band  widths  are  determined  by  the  adjacent  midband  values. 

37.  The  fifth  record  contains  the  list  of  station  numbers  with  I,J 
locations  of  the  stations  in  the  grid  where  model  output  will  be  saved.  In 
this  case,  only  one  station  (MTSA=1)  is  specified  at  location  1=6,  J=7 . 

38.  The  sixth  record  contains  optional  beginning  and  ending  dates  for 
model  execution  and  saving  of  boundary  conditions.  In  this  example,  values 
are  set  at  00000000  and  99999999,  so  the  beginning  and  ending  dates  are  not 
determined  by  these  values,  but  by  the  header  records  on  the  first  and  last 
input  wind  fields,  respectively.  If,  for  example,  one  had  a  year  of  wind 
fields  for  1991  in  a  file  and  only  wanted  to  model  the  month  of  March,  the 
dates  91030100  and  91033121  could  be  specified  for  the  first  and  last  dates, 
respectively.  The  model  would  read  the  wind  input  file  and  begin  when  it 
reached  the  first  date  and  quit  when  it  reached  the  last. 
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39.  The  next  group  of  data  is  the  depths  at  each  grid  point  if  the 
parameter  ISTDEP^O.  If  ISTDEP=1,  deep  water  (999  m)  is  assumed  at  every  grid 
point  and  no  depths  are  read  in.  Depths  are  read  in  for  all  columns  from  the 
top  row  to  the  bottom  row.  The  READ  statement  is  in  Subroutine  READOPT  and 
uses  FORMAT  statement  401.  The  format  of  the  depths  in  the  input  options  file 
must  conform  to  the  format  specified  by  statement  401.  In  this  example,  it  is 
18F4.0. 

40.  The  next  group  of  data  is  the  land-water  boundary  matrix.  This  set 
of  data  is  always  needed  to  specify  which  grid  points  are  land  and  which  are 
water.  It  is  also  used  to  designate  grid  points  at  which  boundary  information 
is  expected  from  a  previous  model  application.  This  example  is  for  an  initial 
application  and  so  no  boundary  points  are  designated  for  input  of  information. 
A  value  of  0  represents  land  and  a  value  of  1  water.  Values  are  read  in  for 
all  columns  from  the  top  row  to  the  bottom  using  FORMAT  statement  957  in 
subroutine  READOPT.  In  this  example,  the  format  is  1811. 

41.  If  the  present  model  application  is  being  used  to  generate  boundary 
information  for  a  smaller  grid  nested  within  the  larger,  the  number  of  points 
defining  the  boundary  of  the  smaller  grid  within  the  larger  is  needed.  In 
this  example,  10  water  points  in  the  present  grid  define  the  water  boundary  of 
a  finer  grid  over  the  southwestern  portion  of  the  lake.  These  points  are 
shown  by  circles  in  Figure  3. 

42.  The  next  two  records  contain  the  column  (I)  and  row  (J)  locations 
of  these  points,  respectively.  The  order  in  which  the  points  are  specified 
can  be  either  clockwise  from  the  western  boundary  of  the  lake  to  the  southern 
boundary,  or  counterclockwise  from  the  southern  boundary  to  the  west  side  of 
the  lake.  Either  way,  the  I  values  must  be  in  the  same  order  as  the  J's,  so 
that  paired  I's  and  J's  define  the  proper  grid  points. 

43.  The  only  other  input  needed  for  this  application  is  the  input  wind 
fields.  An  example  of  the  winds  for  one  time  period  is  shown  below.  The 
first  record  is  the  date/time  group  given  as  four  groups  of  two  digits 
defining  the  year,  month,  day,  and  hour  (from  00  to  23).  The  next  group  of 
records  contain  the  wind  speed  at  all  grid  points  ordered  the  same  as  the 
depths  and  land-water  matrix;  that  is,  values  *"01'  all  columns  starting  with 
the  top  row  and  ending  with  the  bottom.  The  wind  directions  follow  the  wind 
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speeds.  Successive  sets  of  wind  fields  are  supplied  every  NHR  hours  during 
the  simulation. 
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General  Guidance  on  Model  Use 


44.  The  application  of  this  model  begins  with  a  definition  of  what  is 
needed  from  the  model  and  why  it  is  being  used.  The  model  can  be  used  over  a 
range  of  applications;  e.g.,  from  hindcasting  wave  information  on  an  oceanic 
scale  over  decades  to  short-term  simulation  of  events  in  a  bay.  Defining  the 
problem  will  help  set  the  regional  coverage  and  other  input  options.  Note 
that  waves  arriving  at  the  coast  may  have  been  generated  by  storms  hundreds  or 
even  thousands  of  miles  away.  If  the  model  coverage  does  not  extend  to  these 
generation  areas,  this  wave  energy  will  not  be  represented.  Thus  the  first 
step  is  to  define  the  boundaries  of  the  region  which  needs  to  be  modeled. 

45.  Once  these  boundaries  are  established,  a  grid  of  latitude  and 
longitude  lines  is  overlaid  to  define  the  points  (crossings  of  lines  of 
latitude  and  longitude)  at  which  calculations  will  be  made.  The  spacing 
between  lines  should  be  small  enough  to  resolve  geographic  or  bathymetric 
features  of  interest  near  land,  but  not  so  small  that  a  large  number  (on  the 
order  of  100)  of  grid  lines  in  each  direction  results.  Very  large  grids 
require  large  amounts  of  computer  memory  and  computational  time.  The  nested 
grid  feature  of  the  model  should  be  used  when  it  is  required  to  cover  a  large 
region  and  to  also  have  fine  resolution  of  land  boundaries  or  bathymetry  near 
land . 

46.  The  time  step  is  the  next  parameter  to  specify  after  the  grid 
spacing  is  determined.  The  time  step  should  generally  be  as  large  as  possible 
without  exceeding  the  time  it  takes  for  the  fastest  traveling  wave  energy  to 
cross  one  grid  cell.  Since  the  fastest  traveling  waves  are  those  of  the 
lowest  frequency,  a  first  estimate  of  the  time  is  obtained  by  dividing  the 
grid  cell  size  by  the  group  speed  of  the  wave  in  the  lowest  frequency  band. 

The  resulting  time  in  seconds  is  a  safe  value  to  use.  Typically  a  value  close 
to  this  which  is  an  even  multiple  of  the  output  time  is  used.  A  common 
symptom  of  choosing  too  large  a  time  step  is  excessive  growth  of  wave  energy 
resulting  in  unrealistic  wave  heights.  This  is  generally  obvious,  but  can  be 
checked  by  using  Figure  3.23  or  3.24  in  the  SPM  (US  Army  Engineer  Waterways 
Experiment  Station  1984,  pages  3-49  or  3-50)  with  an  estimate  of  wind  speed 
and  fetch. 
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47 .  The  mid-band  values  of  frequency  for  the  frequency  bands  should  be 
chosen  to  represent  all  periods  of  expected  wave  energy.  If,  for  example,  the 
problem  concerns  wave  energy  on  a  small  lake,  long-period  (  >  12  sec)  waves 
will  probably  not  be  present;  thus,  there  would  be  no  need  to  have  bands  in 
this  low-frequency  part  of  the  spectrum.  However,  it  is  desirable  to  have  2-3 
bands  lower  than  the  lowest  band  in  which  the  peak  energy  is  expected  in  order 
to  define  the  shape  of  the  spectrum  in  the  low  frequency  or  forward  face 
region . 

48.  The  program  has  a  set  format  for  reading  the  required  input. 
Generally  there  is  no  reason  to  change  this  for  the  options  file.  Sometimes 
it  is  more  convenient  to  change  the  format  for  input  of  wind  information  than 
to  reformat  input  wind  fields.  If  formats  are  changed  or  data  reformatted  for 
input,  the  user  needs  to  be  sure  that  the  input  data  are  consistent  with  the 
read  statements  and  formats  in  the  program.  Input  read  errors  or 
misinterpretation  of  input  parameters  are  generally  the  result  of 
inconsistency  between  input  data  and  read  statements  or  formats. 

49.  In  a  nested  grid  application,  it  is  assumed  that  the  time  interval 
between  input  wind  fields  for  the  sub-grid  is  the  same  as  that  at  which 
boundary  information  was  saved  in  the  large  grid  application.  This  is  checked 
in  the  program  by  comparing  dates  in  the  header  records  of  wind  fields  and 
boundary  information,  and  will  stop  execution  if  the  dates  don't  match. 
Consistency  of  this  input  should  be  checked  prior  to  application  to  avoid  this 
error.  Also  note  that  the  time  step  will  change  in  a  nested  grid  application 
since  the  grid  spacing  will  be  reduced  for  the  same  climate  of  waves. 

50.  When  forming  the  land/sea  matrix,  it  is  required  that  all  boundary 
points  on  the  circumference  of  the  grid  be  zeros.  This  is  naturally  the  case 
when  a  water  body  is  surrounded  by  land.  If  land  does  not  form  a  natural 
boundary,  the  grid  boundary  should  be  chosen  to  realistically  represent  the 
fetch  length  to  the  region  of  interest.  Likewise,  if  wave  generation  and 
propagation  from  a  water  region  will  not  affect  the  region  of  interest,  it  can 
be  considered  land  and  thus  reduce  the  number  of  calculations. 

51.  Parameter  statements  are  used  in  the  FORTRAN  Code  to  set  the  limits 
on  the  dimensioned  arrays  in  the  program.  These  values  need  to  be  changed  to 
be  appropriate  for  each  application.  The  values  set  are  the  maximum  number  of 
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grid  columns,  grid  rows,  frequency  bands,  angle  bands,  output  boundary  points, 
and  input  boundary  points. 

52.  After  an  application  is  completed,  the  results  should  be  examined 
to  see  if  they  are  reasonable  as  judged  by  nearby  measurements  or 
approximations  for  typical  wind  speeds  and  fetches.  Generally,  a  numerical 
model  is  applied  as  an  interpolation  tool  to  provide  information  in  space  and 
time  where  it  is  lacking.  An  attempt  should  be  made  to  verify  the  model 
results  against  other  available  information  in  order  to  validate  the  results 
and  provide  some  indication  of  accuracy.  Assistance  in  application  of  the 
model  and  interpretation  of  results  is  available  from  the  WIS  project  office. 
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APPENDIX  A:  GENERAL  DESCRIPTION  OF  SUBROUTINES  INITSPR  AND  SPROP 


1.  The  upgrade  of  WISWAVE  I.O  to  include  the  effects  of  refraction  and 
shoaling  in  shallow  water  was  accomplished  by  adding  two  subroutines.  The 
first  is  subroutine  INITSPR,  which  calculates  variables  dependent  on  depth 
that  are  used  in  subroutine  SPROP  to  propagate  wave  energy  from  point  to 
point.  Subroutine  INITSPR  is  called  from  the  main  part  of  the  program  when 
depths  are  defined  initially,  or  when  subsequently  updated  as  in  a  storm  surge 
or  large  tidal  amplitude  case. 

2.  Calculations  in  INITSPR  are  done  within  four  nested  DO  LOOPS.  The 
outermost  loop  is  related  to  the  frequencies  in  the  discrete  spectrum 
(k=l,nfrq),  the  two  intermediate  loops  give  position  (i=2,nx-l,  j=2,ny-l),  and 
the  inner  loop  is  related  to  the  angle  bands  in  the  discrete  spectrum 
(l=l,nang).  Notation  of  variables  and  symbols  below  is  the  same  as  that  used 
in  the  subroutines.  Values  are  calculated  for  five  variables,  each  a  function 
of  frequency,  position,  and  angle. 

3.  The  concept  of  "backtracking"  wave  rays  is  used  to  calculate  these 
variables.  For  a  given  frequency  and  position  (value  of  k,i,j),  a  wave  ray, 
with  a  specified  angle,  can  be  backtracked  from  the  point  using  the  ray 
trajectory  equation 


d  e 

d  s 


1 

ca{i,  j ,  k) 


(i,j)  sin(ll) 


(i,j)  cos(JJ)] 
dy 


(Al) 


for  one  time  step.  The  phase  velocity  is  represented  by  "c,"  and  "II" 
references  the  backtracking  ray.  The  values  dc/dx  and  dc/dy  are  calculated 
using  the  average  value  of  c  over  two  grid  cells  in  the  x  and  y  directions 
respectively.  This  ray  path  is  determined  by  the  initial  angle,  with  respect 
to  the  coordinate  system ,  of  a  tangent  to  the  ray  at  the  (i,j)  point,  and  the 
slope  of  the  bathymetry  defined  by  the  differential  in  depth  of  points  on 
either  side  of  (i,j)  in  the  north-south  and  east-west  directions.  The  arc 
length  traced  by  the  ray  is  part  of  a  circular  path. 
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4.  For  a  given  frequency  and  position,  the  ray  arrival  angle  is  varied 
through  all  direction  bands  relative  to  the  bottom  slope  angle  (angsl)  to 
allow  computation  of  the  change  in  angle  (dth) ,  component  displacements  (dpx, 
dpy)  in  the  x,  y  directions,  and  combined  component  refraction,  shoaling 
coefficients  (ccgrjx,  ccgrjy).  These  are  the  variables  transferred  to 
subroutine  SPROP. 

5.  Subroutine  SPROP  examines  the  amount  of  energy,  in  x,  y  component 
form,  leaving  the  point  (i,j)  and  arriving  at  (i,j)  from  the  upstream  grid 
point  accounting  for  refraction  and  shoaling.  A  more  detailed  description  of 
the  calculations  is  provided  by  comments  in  the  coding  of  these  two 
subroutines  (see  Appendix  B) . 
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APPENDIX  B;  SUBROUTINES  ADDED  IN  UPGRADE  TO  WISWAVE  2.0 


subroutine  initspr 

c  This  subroutine  is  called  once  from  the  main  program  if  the  shallow  water 
c  option,  IDRP~0  is  input.  It  is  not  called  if  depths  are  set  to  999  m  by 
c  setting  ISTDEP=1,  or  if  depths  are  input  but  IDPRP=1 .  If  depths  are  updated 
c  to  reflect  changing  water  levels,  it  is  called  from  the  main  program  each 
c  time  an  update  occurs . 

include  'paral.inc' 

c  The  include  statement  reads  a  file  (paral.inc)  which  is  a  fortran  parameter 
c  statement  used  to  set  the  values  of  idmn, jdmn, if , ia  for  dimensioning 
c  arrays. 

common  /smod/  dpx(idmn, jdmn, if , ia) ,  dpy ( idmn , j dmn , if , ia) , 

2dth( idmn , j  dmn , if , ia) , ccgrj  x( idmn , j  dmn , if , ia) , ccgrjy ( idmn , 

3j  dmn , if , ia) 

c  dpx,  dpy  are  component  distances,  normalized  by  the  grid  cell  size, 
c  traveled  by  a  backtracked  ray  in  one  time  step. 

c  dth  is  the  change  in  refraction  angle  from  upstream  to  point  (i,j). 
c  ccgrjx,  ccgrjy  are  related  to  refraction,  shoaling  coefficients. 

common  /ccg/  ca( idmn , j dmn , i f) , cga ( idmn , jdmn , if ), delx(j dmn) 
c  ca,  cga  are  respectively  the  phase  and  group  velocities  at  grid  points, 
c  delx  is  latitude  dependent  due  to  the  convergence  of  meridians. 

common  /a3/  g , pi , dl , nang , nf rq , nx , ny , twopi , dt , raddeg 
c  g=9 . 8  m/sec*'*2  ,  pi=3 . 19159  ,  dl=distance  between  grid  points  (km),  nang“#  of 
c  angle  bands,  nfrq=//  of  frequency  bands,  nx=#  of  columns  ingrid,  ny=//  of 
c  rows  in  grid,  twopi=2*pi,  dt=time  step  (sec),  raddeg=360/twopi , 

dimension  dcgdx( idmn, j  dmn) , dcgdy( idmn, j  dmn) , dcdx( idmn, j  dmn) , 

2dcdy ( idmn , j  dmn) 

c  dcgdx,  dcgdy  are  respectively  the  change  in  group  velocity  in  the  x,y 

c  directions  about  point  (i,j)  using  centered  differences.  dcdx,dcdy  are 

c  equivalent  values  for  the  phase  velocity. 

common  /o5/  ibound( idmn , j dmn) 

c  ibound  is  the  land- water  boundary  matrix,  which  can  include  boundary 
c  extrapolation,  interpolation,  or  coincident  point  (same  point  in  parent  and 
c  nested  grids)  indicators. 

common  /angs/  ang( ia) , sina( ia) , cosa( ia) , islp ( idmn , j dmn) 
c  ang(l)  is  the  center  angle  in  radians  of  angle  band  1,  sina  and  cosa  are 
c  respectively  the  sine  and  cosine  of  ang(l),  islp(i,j)  is  an  index  (either  0 

c  or  1)  indicating  which  way  the  bottom  slopes  at  i,j,  more  details  below. 

common  /depth/  dep ( idmn , j dmn) 
c  The  values  of  depth  (m)  at  each  grid  point. 

common  /band/  iua( ia) , juaa( ia) 
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c  Indices  (0  or  +,-  1)  used  to  locate  the  appropriate  upstream  grid  point 
c  from  i,j  for  each  angle  band. 


ep=0. 0000001 

c  A  small  value  added  to  the  denominator  to  avoid  dividing  by  zero, 
ainc-twopi/nang 

c  The  increment  between  angle  bands  in  radians 
alim=2*ainc 

c  A  limit  on  the  change  in  angle  due  to  refraction  between  an  upstream  point 
c  and  i , j . 

do  80  l=l,nang 
ang (l)=(l-l)*ainc 
sina(l)=sin(ang(l) ) 
cosa(l)=cos(ang(l)) 

80  continue 

c  Calculate  angles  (radians)  and  sines,  cosines  for  all  angle  bands 

do  86  i=l,nang 

x=cosa(i) 

y=sina(i) 

iua(i)=-sign(1.0,x) 
juaa( i)=-sign(1.0,y) 
if  (abs (x) . It . 0 . 0001 )  iua(i)=0.0 
if  (abs (y) . 1 t . 0 . 0001 )  juaa(i)=0.0 
86  continue 

c  calculate  the  closest  upstream  i,j  point  from  (i.j)  for  each  angle  band. 
curvlim=l . /( 2 . *dl ) 

c  This  is  a  limit  on  the  arc  length  generated  by  the  change  in  refraction 
c  angle.  It  is  equal  to  1/radius  of  curvature  which  is  taken  as  2*dl . 

do  50  k=l,nfrq 
do  60  i=l , nx 
do  60  j  =1 , ny 
do  60  1=1 , nang 
dpx( i , j , k , 1 ) =0 . 
dpy ( i , j , k , 1 ) =0 . 
dth( i , j , k , 1 )=0 . 
ccgrj  x(i,j ,k,l)=l. 
ccgrjy(i,j ,k,l)=l. 

60  continue 

c  Initialize  values  of  arrays 

do  1  i=2,nx-l 
do  1  j=2  ,  ny-  1 

dcgdx( i,j)=(cga(i+l,j,k)-cga(i-l,j,k))/(2. *delx( j ) ) 
dcgdy( i ,j)  =  (cga(i,j+l,k)- cga( i , j - 1 ,k) )/(2 .*dl) 
dcdx(i,j)=(ca(i+l,j,k)-ca(i-l,j,k))/(2. *de lx( j ) ) 
dcdy(i,j)=(ca(L,j+l,k)-ca(i,j-l,k))/(2.*dl) 
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1  continue 

c  Calculate  differentials  of  group  and  phase  velocity. 

do  2  i-=2,nx-l  ; 

do  2  j “2,07-1  ■ 

if  (ibound(i , j ) .ne . 1)  go  to  2 

c  Skip  to  end  of  loop  if  not  a  water  point.  Depths  are  set  to  0  for  land 
c  points  in  subroutine  rdopt 

dhx-dep ( i - 1 , j ) - dep ( i+1 , j ) 
dhy-dep ( i , j - 1 ) - dep ( i , j +1 ) 
angsl=atan2 (dhy , dhx+0 . 0000001) 

c  angsl  is  the  slope  angle  (positive  toward  shallow  water) 

islp(i , j )=0 

if  (abs (cos (angsl )). gt . abs (sin(angsl) ) )  islp(i,j)=l 
c  islp=l  for  angles  >  315  and  <  45,  and  >  135  and  <  225  deg  cartesian  system, 
c  otherwise,  islp=0. 

do  3  l=l,nang 
ll=l+nang/2 

if  (ll.gt.nang)  ll=ll-nang 

c  1  and  11  are  indices  for  opposite  angle  bands,  11  references  the 
c  "backtracking"  ray. 

cgav“cga( i , j , k)+  0 . 5*dcgdx( i , j )*cga(i , j , k)*cosa( ll)*dt 
2  +  0 . 5*dcgdy ( i , j )*cga( i , j , k)*sina( ll)*dt 
c  cgav  is  the  average  group  velocity  calculated  from  points  around  i,j. 

dels“Cgav*dt 

c  distance  energy  will  travel  in  one  time  step  at  the  average  group  speed. 

dthds“(dcdx(  i  ,j)*sina(ll)-dcdy(i,j)’*^cosa(ll))/(ca(i,j  ,k)+ep) 
c  This  is  the  ray  trajectory  equation.  It  calculates  the  change  in  angle  due 
c  to  refraction  along  an  arc  length. 

if  (abs (dthds) . gt . curvlim)  dthds=curvlim*dthds/abs (dthds) 
c  Limit  dthds  to  curvlim,  retaining  the  appropriate  sign  (+  or  -). 

thtav=ang( 11 )+0 . 5*dthds*dels 

c  An  average  angle,  based  on  initial  value  +  1/2  the  change 

phi=1.0/(sqrt(1.0-(0. 5*dels*dthds)**2) ) 
c  1/ratio  of  arc  length  to  chord  length.  Chord  length  approximated  with  a 
c  Taylor  expansion  of  2r  sin  (1/2  theta) 

dpx( i , j , k , 1 )=dels*cos ( thtav) /phi 
dpy ( i , j , k , 1 ) =de ls*s in ( thtav) /phi 
dpx( i , j , k , 1 )=abs (dpx ( i , j , k , 1) )/delx( J ) 
dpy ( i , j , k , 1 ) =abs ( dpy ( i , j , k . 1 ) ) /dl 

c  Calculate  fraction  of  cell  size  in  x,y  direction  traveled  by  ray  in  one 
c  time  step 
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ii“i+iua(l) 
cgupx=cga( ii , j ,k) 
cgrx-cgupx/cga(i , j ,k) 
j j=j +juaa(l) 
cgupy=cga(i , j j ,k) 
cgry=cgupy/cga(i , j ,k) 

c  Ratio  of  component  upstream  group  speeds  to  group  speed  at  (i,j) 
c  Solve  for  intersection  with  boundary 

c  Following  section  considers  rays  backtracked  from  (i,j)  at  angles  of  each 
c  angle  band  and  determines  the  originating  angle  based  on  the  bottom  slope 
c  and  phase  speeds  using  Snell's  Law.  Then  the  change  in  angle,  dth ,  is 
c  calculated  along  with  the  proper  ratios  of  group  speed  for  use  in 
c  subroutine  sprop . 

cl=ca(i , j ,k) 
c2=ca(ii, j ,k) 
rx=cl/c2 
c3=ca(i ,j j ,k) 
ry=cl/c3 

c  Ratio  of  phase  speed  at  (i,j)  to  components  of  phase  speed  upstream. 
angz=ang ( 1 ) - angs 1 

c  Difference  between  band  angle  and  slope  angle. 

if  (ang2.gt.pi)  angz=angz - twopi 
if  (angz  .  It . -pi)  angz=angz-H twopi 
c  Express  angle  as  value  between  +  or  -  180  deg. 

if  (abs (angz ). gt . pi/2 . )  rx=l./rx 
if  (abs  (angz)  .  gt .  pi/2  .  )  ry"=l./ry 
c  Invert  ratio  of  phase  speeds  if  angle  >  90  deg. 

if  (angz . gt . pi/2 . )  angz^pi  -  angz 
if  (angz . It . -pi/2 . )  angz=-pi  -  angz 
c  Express  angle  as  value  between  -t-  or  -  90  deg. 

sinz=sin(angz) 
if  (rx.gt.2.0)  rx=2.0 
c  Limit  max.  ratio  of  phase  speeds  to  2 

if  ( rx . gt . 1 . )  go  to  76 
c  Two  options  for  x  component. 

c  If  speed  increasing  as  approaching  (i,j)  go  to  76, 
c  otherwise  speed  decreasing,  do  section  below. 

angO=asin(rx  *sinz) 

c  Snell's  Law  s in( shal low) /s in(deep)=c ( shallow) /c ( deep) . 

bx=abs (cos ( angz) / (cos ( ang0)+ep) ) 
c  Ratio  of  ray  separation  in  deep  to  shallow  water, 
c  bx  and  by  are  related  to  refraction  coefficient 


if (abs(sinz) . gt . 0 . 99985)  bx=l . 0 
c  Set  bx  for  90  or  270  deg  incidence. 

if  ( islp( i , j ) . eq . 1)  dth(i , j ,k, l)=angO-angz 
c  Calculate  change  in  angle. 

go  to  61 

c  Do  y  component 

c  Arrive  here  if  speed  increasing  as  approaching  (i,j) 

76  angO“asin(sinz/rx) 

bx^abs (cos (ang0)/(cos (angz )+ep) ) 

if  ( islp( i , j ) . eq . 1)  dth(i , j ,k,l)=angz-angO 

c  Two  options  for  y  component 

61  if  (ry.gt.2.0)  ry=2.0 
if  (ry.gt.l.)  go  to  62 
angO=asin(ry*sinz) 

by=abs (cos (rngz) /(cos (ang0)+ep) ) 
if (abs ( s inz ). gt . 0 . 99985 )  by=1.0 
c  Set  by  for  90  or  270  deg  incidence 

if  ( islp ( i , j ) . ne . 1 )  dth( i , j , k , l)=angO-angz 
go  to  63 

c  Arrive  here  if  speed  increasing  as  approaching  (i,j) 

62  angO”asin(sinz/ry) 
by=abs(cos(angO) /(cos(angz)+ep) ) 
if (abs (sinz) . gt . 0 . 99985)  by=1.0 

c  Set  by  for  90  or  270  deg  incidence 

if  ( islp ( i , j ) . ne . 1 )  dth( i , j , k , l)=angz-ang0 

63  if  (dth( i , j , k , 1 ) . gt . alim)  dth( i , j , k , l)=alim 
if  (dth( i , j , k , 1) . It . - alim)  dth( i , j ,k , 1)=- alim 

c  Limit  change  in  angle  to  two  angle  increments. 

if  (bx.gt.l.)  bx  =  l./cgrx 
if  (by.gt.l.)  by  =  l./cgry 
ccgrj x( i , j , k , 1 )=cgrx  *bx 
ccgrjy(i , j ,k, l)=cgry  *by 

c  ccgrjx,  ccgrjy  are  related  to  combined  shoaling,  refraction  coefficient. 

3  '.;ontinue 

c  End  of  angle  loop 
2  continue 

c  End  of  space  loop 

50  continue 

c  End  of  frequency  loop 

c  The  following  variables  are  calculated  above  and  used  in  subroutine  sprop 
c  Variable  Function  of 

c- . - . . . . 
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c  ispl  i,j 

c  dpx ,  dpy  i , j , k , 1 

c  ccgrjx,  ccgrjy  i,j,k,l 

c  dth  i , j , k , 1 

return 

end 

subroutine  sprop(k) 

c  This  subroutine  is  called  from  the  main  program  within  loop  2000,  which  is 
c  a  loop  on  frequency.  It,  in  turn,  is  in  loop  1001,  which  is  a  loop  on  time 
c  steps.  Loops  in  this  subroutine  are  on  space  and  angles,  thus  computations 
c  in  this  subroutine  are  done  for  every  time  step,  frequency,  grid  point  and 
c  angle.  The  purpose  of  the  subroutine  is  to  propagate  energy  in  shallow 
c  water  accounting  for  refraction  and  shoaling.  The  equivalent  subroutine 
c  for  deep  water  is  dprop . 

include  'paral.inc' 

c  The  include  statement  reads  a  file  (paral.inc)  which  is  a  fortran  parameter 
c  statement  used  to  set  the  values  of  idmn , j dmn , if , ia  for  dimensioning 
c  arrays. 

common  /smod/  dpx(idmn, jdmn, if , ia) ,  dpy ( idmn , j dmn , if , ia) , 

2  dth( idmn , j  dmn , if , ia) , ccgrjx( idmn, j  dmn, if , ia) , ccgrjy ( idmn , 

3  j  dmn, if , ia) 

c  dpx,  dpy  are  component  distances,  normalized  by  the  grid  cell  size, 
c  dth  is  the  change  in  refraction  angle  from  upstream  to  point  (i,j). 
c  ccgrjx,  ccgrjy  are  related  to  refraction,  shoaling  coefficients. 

common  /ccg/  ca ( idmn , j dmn , if ) , cga( idmn, j dmn, if ) , delx( j dmn) 
c  ca,  cga  are  respectively  the  phase  and  group  velocities  at  grid  points. 

common  /a3/  g , pi ,xn, dl ,nang,nfrq ,nx ,ny , twopi , fac , dt , raddeg 
c  g=9 . 8  m/sec**2 ,  pi=3 . 14159 , dl=distance  between  grid  points  (km),  nang=//  of 
c  angle  bands,  nfrq=y/  of  frequency  bands,  nx=//  of  columns  ingrid,  ny=(/  of 
c  rows  in  grid,  twopi=2*pi,  dt=time  step  (sec),  raddeg=360/twopi . 

dimension  eoutx( idmn, jdmn, ia) , eouty( idmn , j dmn , ia) , 
leinx ( idmn , j  dmn , ia) , einy ( idmn , j  dmn, ia) , e inf x( idmn , j  dmn , ia) , 

2e infy ( idmn , j  dmn , ia) 

c  eoutx.eouty  are  components  of  energy  leaving  point  (i,j)  at  angle  1. 
c  einx.einy  are  components  of  energy  entering  point  (i,j)  at  angle  1. 
c  einfx,einfy  are  components  of  energy  for  all  angles. 

common  /angs/  ang( ia) , s ina ( ia) , cosa( ia) , islp ( idmn , j dmn) 
c  ang(l)  is  the  center  angle  in  radians  of  angle  band  1,  sina  and  cosa  are 
c  respectively  the  sine  and  cosine  of  ang(l),  islp(i,j)  is  an  index  (either  0 
c  or  1)  indicating  which  way  the  bottom  slopes  at  i,j,  more  details  below. 

common  /cen/  e ( idmn , j dmn , if , ia) , enxt (nbn , if , ia) , elst (nbn , if , ia ) 
c  e  is  the  frequency  spectrum  at  (i,j)  after  propagation.  enxt  and  elst  are 
c  not  used  in  this  subroutine,  see  main  program. 


B6 


conunon  /o5/  ibound(idmn, jdmn) 
c  ibound  is  the  land- water  boundary  matrix. 

common  /band/  iua( ia) , j uaa( ia) 

c  Indices  (0  or  +,-  1)  used  to  locate  the  appropriate  upstream  grid  point 
c  from  i , j . 

ainc=twopi/nang 

c  The  increment  between  angle  bands  in  radians 
c  Entering  x- component  propagation  section 

c  Determine  x  component  of  energy  leaving  (i,j)  and  arriving  at  (i,j)  from 
c  point  upstream  for  all  angle  bands  at  one  frequency 
do  1  i=l,nx 
do  1  j  =  1 ,  ny 

if  (ibound(i , j ) .ne . 1)  go  to  1 
do  2  l=l,nang 

eoutx(i,j ,l)=e(i,j ,k, 1 )*dpx( i , j ,k, 1) 
ii=i+iua(l) 

e inx (i,j ,l)=e(ii,j ,k,l) 

2  continue 
1  continue 

c  Entering  theta-propagation,  x  component 
do  5  i=2,nx-l 
do  5  j=2,ny-l 

if  ( ibound( i , j ) . ne . 1)  go  to  5 
c  If  not  a  water  point  go  to  next  point 

do  6  l=I,nang 
einfx(i, j , 1)=0. 

6  continue 

c  Initialize  energy  summation  array  at  all  angles  to  zero 
do  7  1=1 , nang 

angz=ang(l)-^dth(i,  j  ,  k ,  l)*islp(  i ,  j  ) 
if (abs (angz) . It . 0 . 02)  angz=0.0 
if (angz . gt . twopi)  angz=angz- twopi 
if (angz . It . 0 . 0)  angz=angz+twopi 
ll=angz/ainc+l 
if(ll.lt.l)  angz=angz+twopi 
if (11. It. 1)  ll=ll+nang 
if ( 11 . gt . nang)  ll=ll-nang 
lh=ll-H 

if ( Ih . gt . nang)  lh=lh-nang 
xl=(angz - (11  - l)*ainc)/ainc 

einfx(i,j ,ll)=einfx(i,j , ll)+ccgrjx( i , j , k , l)*einx( i ,j ,1)*(1. -xl) 

2  *  dpx(i, j ,k,ll) 

einfx( i , j , lh)=einfx( i , j , lh)+ccgrjx( i , j , k , l)*einx( i , j , l)*xl 
1  *  dpx( i , j ,k, Ih) 

/  continue 
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c  Determine  the  amount  of  energy  arriving  at  (i,j)  from  adjacent  upstream 
c  angle  bands . 

do  8  1=1 , nang 

e(i,j ,k,l)=e(i,j ,k,l) -eoutx(i , j , l)+einfx(i , j , 1) 
if  (e(ij  ,k,l)  .It.O.  )  e(i,j  ,k,l)=0. 

8  continue 

c  Calculate  new  energy  after  propagation  in  each  angle  band  at  frequency  k 
c  and  point  ( i , j ) 

5  continue 

c  Completion  of  space  loop 

c  Entering  y-component  propagation  section 

c  Determine  y  component  of  energy  leaving  (i,j)  and  arriving  at  (i,j)  from 
c  point  upstream  for  all  angle  bands  at  one  frequency 

do  17  i=l,nx 
do  17  j=l,ny 

if  ( ibound( i , j ) . ne . 1)  go  to  17 
c  If  not  a  water  point  go  to  next  point 

do  9  1=1 , nang 

eouty(i, j ,l)=e(i,j , k , 1 )*dpy ( i , j ,k,l) 

j j=i+juaa(l) 

einy(i,j ,l)=e(i,jj ,k,l) 

9  continue 
17  continue 

c  Entering  theta-propagation,  y  component 
do  3  i=2,nx-l 
do  3  j=2,ny-l 

if  ( ibound(i , j ) , ne . 1)  go  to  3 
do  16  l=l,naiig 
einfy ( i , j , 1)=0 . 

16  continue 

do  27  1=1, nang 

angz=ang ( 1 ) +d th ( i , j , k , 1 ) * ( 1 -  is Ip ( i , j ) ) 

if (abs (angz) . It . 0 . Oz )  angz=0 . 

i f ( angz . gt . twop i )  angz=angz - twop i 

if (angz . It . 0 . )  angz=angz+twopi 

ll=angz/ainc+l 

if(ll.lt.l)  ll=ll+nang 

if ( 11 . gt . nang)  ll=ll-nang 

lh=ll-(-l 

if ( Ih . gt . nang)  lh=lh-nang 
xl=(angz - ( 11  - 1 )*ainc)/ainc 

einfy(i,j , ll)=einfy( i , j , ll)+ccgr jy( i , j ,k, l)*einy(i , j ,1)*(1. -xl) 
2*dpy(i,j ,k,ll) 

einfy (i,j  ,lh)=einfy(i,j  ,lh)+ccgrjy(i,j  ,k,l)*einy(i,j  ,l)''^xl 
1'^  dpy  (  i  ,  j  ,  k  ,  Ih) 
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27  continue 

c  Determine  the  amount  of  energy  arriving  at  (i,j)  from  adjacent  upstream 
c  angle  bands . 

do  18  l=l,nang 

e(i,j ,k,l)=e(i,j ,k,l)- eouty ( i , j , l)+einfy ( i , j , 1) 
if  (e(i,j ,k,l) .It.O.)  e(i,j.k,l)=0. 

18  continue 

c  Calculate  new  energy  after  propagation  in  each  angle  band  at  frequency  k 
c  and  point  ( i , j ) 

3  continue 

c  Completion  of  space  loop 

return 

end 
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